import math

def spherical_interpolation(lat1, lon1, lat2, lon2, num_points):
    """
    在两个航点之间进行球面线性插值
    :param lat1: 起始点纬度(度)
    :param lon1: 起始点经度(度)
    :param lat2: 终点纬度(度)
    :param lon2: 终点经度(度)
    :param num_points: 返回的总点数(包括起始点和终点)
    :return: 插值点列表[(lat, lon), ...]
    """
    if num_points < 2:
        return [(lat1, lon1), (lat2, lon2)][:num_points]
    
    # 角度转弧度
    lat1_rad = math.radians(lat1)
    lon1_rad = math.radians(lon1)
    lat2_rad = math.radians(lat2)
    lon2_rad = math.radians(lon2)
    
    # 转换为三维笛卡尔坐标（单位球面）
    x0 = math.cos(lat1_rad) * math.cos(lon1_rad)
    y0 = math.cos(lat1_rad) * math.sin(lon1_rad)
    z0 = math.sin(lat1_rad)
    
    x1 = math.cos(lat2_rad) * math.cos(lon2_rad)
    y1 = math.cos(lat2_rad) * math.sin(lon2_rad)
    z1 = math.sin(lat2_rad)
    
    # 计算向量夹角
    dot = x0*x1 + y0*y1 + z0*z1
    dot = max(min(dot, 1.0), -1.0)  # 防止浮点误差
    theta = math.acos(dot)
    
    # 处理重合点
    if theta < 1e-10:
        return [(lat1, lon1)] * num_points
    
    results = []
    for i in range(num_points):
        t = i / (num_points - 1)
        
        # 球面线性插值公式
        a = math.sin((1 - t) * theta) / math.sin(theta)
        b = math.sin(t * theta) / math.sin(theta)
        
        x = a * x0 + b * x1
        y = a * y0 + b * y1
        z = a * z0 + b * z1
        
        # 转换回经纬度
        lat_rad = math.atan2(z, math.sqrt(x**2 + y**2))
        lon_rad = math.atan2(y, x)
        
        # 弧度转角度
        lat_deg = math.degrees(lat_rad)
        lon_deg = math.degrees(lon_rad)
        
        # 规范化经度到[-180, 180]范围
        lon_deg = (lon_deg + 180) % 360 - 180
        
        results.append((lat_deg, lon_deg))
    
    return results

# 示例使用
if __name__ == "__main__":
    # 北京首都机场 (40.0800°N, 116.5900°E)
    # 上海浦东机场 (31.1433°N, 121.8050°E)
    points = spherical_interpolation(40.0800, 116.5900, 31.1433, 121.8050, 5)
    
    print("inserted points:")
    for i, (lat, lon) in enumerate(points):
        print(f"pos{i+1}: lat={lat:.6f} deg, lng={lon:.6f} deg")